#---------------------------------------------------------
# Are Intermediary Constraints Priced?
# Du, Hébert, Huber
# The Review of Financial Studies 2022
#---------------------------------------------------------

#---------------------------------------------------------
# Calculate summary statistics and Sharpe ratios in tables
# To be used in final formatted output
#---------------------------------------------------------

#---------------------------------------------------------
# Set-up
#---------------------------------------------------------

if (exists('script_path')) {
  rm(list = setdiff(ls(), 'script_path'))
} else {
  args = commandArgs(trailingOnly = T)
  if (length(args) > 0) {
    script_path <- args[1]  
    rm(list = setdiff(ls(), 'script_path'))
  } else {
    rm(list = ls())
    script_path <- '~/Dropbox/ForwardArbitrage/CodeReplication/'
  } 
}

# Inputs
source(paste0(script_path, 'RCode/ArbFunc.R'))
load(paste0(script_path, 'OutputInterim/Basis_and_Returns.RData')) 
load(paste0(script_path, 'OutputInterim/All_Rates_Dates_Basis_Returns.RData')) 

# Initialization
currency_names6 <- generate_currency_names(number_of_currencies = 6)
basis_matrix_row <- c('sample_size', 'mean', 'sd_mean')
return_noSE_row <- c('sample_size', 'mean', 'sd_mean', 'Sharpe')
return_matrix_row <- c('sample_size', 'mean', 'se_mean', 't_mean', 'sd_mean', 'Sharpe', 'se_Sharpe', 't_Sharpe')

# Create new summary timeline
basis_and_returns[, covid_period := ifelse(Date >= covid_start_date, covid_order[2],
                                          ifelse(Date >= post_crisis_date, covid_order[1], NA))]
basis_and_returns[, slr_period := ifelse(Date >= slr_start_date, slr_order[2],
                                           ifelse(Date >= post_crisis_date, slr_order[1], NA))]

#---------------------------------------------------------
# Single currency returns - top 10 currency pairs AND single against USD
# ** Note: annualized return & SR **
# ** Error messages about GMM arise because COVID and SLR sub samples contain NA; can safely ignore **
# ** Warning messages about xtable arise because of float in LaTeX; can safely ignore **
#---------------------------------------------------------

list_single_return <- list()
list_single_mean_star <- list()
list_single_Sharpe_star <- list()

list_single_return_covid <- list()
list_single_mean_star_covid <- list()
list_single_Sharpe_star_covid <- list()

list_single_return_slr <- list()
list_single_mean_star_slr <- list()
list_single_Sharpe_star_slr <- list()

counter <- 1

for (rate in c('ois', 'ibor')) {
  if (rate == 'ois') {
    currency_names <- c(paste(ois_table$fund_currency, ois_table$invest_currency, sep = '_'),
                        paste(single_table$fund_currency, single_table$invest_currency, sep = '_'))
    return_name_list <- c('1M_fwd_1M', '1M_fwd_3M', '3M_fwd_3M')
  } else if (rate == 'ibor') {
    currency_names <- c(paste(ibor_table$fund_currency, ibor_table$invest_currency, sep = '_'),
                        paste(single_table$fund_currency, single_table$invest_currency, sep = '_'))
    return_name_list <- c('1M_fwd_3M', '3M_fwd_3M')
  }
  
  for (return_tenor in return_name_list) {
    log_return_names <- paste(currency_names, rate, 'log_return', return_tenor, sep = '_')

    # regular 3-period
    temp_matrix <- data.frame(matrix(, nrow = length(return_matrix_row) * length(log_return_names), ncol = length(period_order)))
    colnames(temp_matrix) <- c(period_order)
    rownames(temp_matrix) <- paste(sapply(substr(log_return_names, 0, 7), function(x) rep(x, length(return_matrix_row))),
                                   return_matrix_row)
    for (i in 1:length(log_return_names)) {
      for (p in 1:length(period_order)) {
        temp_matrix[grepl(substr(log_return_names[i], 0, 7), rownames(temp_matrix)), period_order[p]] <- 
          unlist(summary_stats_returns(temp_data_table = basis_and_returns[period == period_order[p]], column_to_summ = log_return_names[i], 
                                       mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                       basis_or_percentage = 'basis'))
      }
    }
    
    temp_title <- paste('Summary stats on single-currency returns of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_single_return[[counter]] <- temp_matrix
    
    star_mean <- matrix(, nrow = length(currency_names), ncol = length(period_order))
    colnames(star_mean) <- period_order
    rownames(star_mean) <- currency_names
    
    star_Sharpe <- matrix(, nrow = length(currency_names), ncol = length(period_order))
    colnames(star_Sharpe) <- period_order
    rownames(star_Sharpe) <- currency_names
    
    # COVID pre/post
    temp_matrix <- data.frame(matrix(, nrow = length(return_matrix_row) * length(log_return_names), ncol = length(covid_order)))
    colnames(temp_matrix) <- c(covid_order)
    rownames(temp_matrix) <- paste(sapply(substr(log_return_names, 0, 7), function(x) rep(x, length(return_matrix_row))),
                                   return_matrix_row)
    for (j in 1:length(log_return_names)) {
      for (q in 1:length(covid_order)) {
        temp_summary <- tryCatch(unlist(summary_stats_returns(temp_data_table = basis_and_returns[covid_period == covid_order[q]], column_to_summ = log_return_names[j], 
                                                         mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                                         basis_or_percentage = 'basis')),
                                 error = function(e) {}) #error from no CHF post-2020
        if (is.null(temp_summary)) {
          temp_matrix[grepl(substr(log_return_names[j], 0, 7), rownames(temp_matrix)), covid_order[q]] <- NA
        } else {
          temp_matrix[grepl(substr(log_return_names[j], 0, 7), rownames(temp_matrix)), covid_order[q]] <- temp_summary
        }
      }
    }

    temp_title <- paste('COVID stats on single-currency returns of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_single_return_covid[[counter]] <- temp_matrix
    
    star_mean_covid <- matrix(, nrow = length(currency_names), ncol = length(covid_order))
    colnames(star_mean_covid) <- covid_order
    rownames(star_mean_covid) <- currency_names
    
    star_Sharpe_covid <- matrix(, nrow = length(currency_names), ncol = length(covid_order))
    colnames(star_Sharpe_covid) <- covid_order
    rownames(star_Sharpe_covid) <- currency_names
    
    # 2015 pre/post
    temp_matrix <- data.frame(matrix(, nrow = length(return_matrix_row) * length(log_return_names), ncol = length(slr_order)))
    colnames(temp_matrix) <- c(slr_order)
    rownames(temp_matrix) <- paste(sapply(substr(log_return_names, 0, 7), function(x) rep(x, length(return_matrix_row))),
                                   return_matrix_row)
    for (j in 1:length(log_return_names)) {
      for (q in 1:length(slr_order)) {
        temp_summary <- try(unlist(summary_stats_returns(temp_data_table = basis_and_returns[slr_period == slr_order[q]], column_to_summ = log_return_names[j], 
                                                         mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                                         basis_or_percentage = 'basis')))
        if (class(temp_summary) == 'try-error') {
          temp_matrix[grepl(substr(log_return_names[j], 0, 7), rownames(temp_matrix)), slr_order[q]] <- NA
        } else {
          temp_matrix[grepl(substr(log_return_names[j], 0, 7), rownames(temp_matrix)), slr_order[q]] <- temp_summary
        }
      }
    }
    
    temp_title <- paste('SLR stats on single-currency returns of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_single_return_slr[[counter]] <- temp_matrix
    
    star_mean_slr <- matrix(, nrow = length(currency_names), ncol = length(slr_order))
    colnames(star_mean_slr) <- slr_order
    rownames(star_mean_slr) <- currency_names
    
    star_Sharpe_slr <- matrix(, nrow = length(currency_names), ncol = length(slr_order))
    colnames(star_Sharpe_slr) <- slr_order
    rownames(star_Sharpe_slr) <- currency_names
    
    for (c in 1:length(currency_names)) {
      star_mean[currency_names[c], ] <- make_stars(values = list_single_return[[counter]][paste(currency_names[c], 't_mean'), ],
                                                   df = list_single_return[[counter]][paste(currency_names[c], 'sample_size'), ])
      star_Sharpe[currency_names[c], ] <- make_stars(values = list_single_return[[counter]][paste(currency_names[c], 't_Sharpe'), ],
                                                     df = list_single_return[[counter]][paste(currency_names[c], 'sample_size'), ])
      star_mean_covid[currency_names[c], ] <- make_stars(values = list_single_return_covid[[counter]][paste(currency_names[c], 't_mean'), ],
                                                   df = list_single_return_covid[[counter]][paste(currency_names[c], 'sample_size'), ])
      star_Sharpe_covid[currency_names[c], ] <- make_stars(values = list_single_return_covid[[counter]][paste(currency_names[c], 't_Sharpe'), ],
                                                     df = list_single_return_covid[[counter]][paste(currency_names[c], 'sample_size'), ])
      star_mean_slr[currency_names[c], ] <- make_stars(values = list_single_return_slr[[counter]][paste(currency_names[c], 't_mean'), ],
                                                         df = list_single_return_slr[[counter]][paste(currency_names[c], 'sample_size'), ])
      star_Sharpe_slr[currency_names[c], ] <- make_stars(values = list_single_return_slr[[counter]][paste(currency_names[c], 't_Sharpe'), ],
                                                           df = list_single_return_slr[[counter]][paste(currency_names[c], 'sample_size'), ])
    }
    
    list_single_mean_star[[counter]] <- star_mean
    list_single_Sharpe_star[[counter]] <- star_Sharpe
    list_single_mean_star_covid[[counter]] <- star_mean_covid
    list_single_Sharpe_star_covid[[counter]] <- star_Sharpe_covid
    list_single_mean_star_slr[[counter]] <- star_mean_slr
    list_single_Sharpe_star_slr[[counter]] <- star_Sharpe_slr
    
    counter <- counter + 1
  }
}
names(list_single_return) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_single_mean_star) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_single_Sharpe_star) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

names(list_single_return_covid) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_single_mean_star_covid) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_single_Sharpe_star_covid) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

names(list_single_return_slr) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_single_mean_star_slr) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_single_Sharpe_star_slr) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

#---------------------------------------------------------
# Average stats for pairs
#---------------------------------------------------------

# Use rates_and_dates because basis_and_returns do not have interest rates
rates_temp <- rates_and_dates[Date >= post_crisis_date & Date <= trading_end_date]

# Average interest rate and basis for the top pairs
list_interest <- list()
counter <- 1
for (rate in c('ois', 'ibor')) {
  if (rate == 'ois') {
    currency_names <- c(paste(ois_table$fund_currency, ois_table$invest_currency, sep = '_'),
                        paste(single_table$fund_currency, single_table$invest_currency, sep = '_'))
    interest_name_list <- c('1M', '3M')
  } else if (rate == 'ibor') {
    currency_names <- c(paste(ibor_table$fund_currency, ibor_table$invest_currency, sep = '_'),
                        paste(single_table$fund_currency, single_table$invest_currency, sep = '_'))
    interest_name_list <- c('3M')
  }
  
  for (interest_name in interest_name_list) {
    temp <- cbind(t(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)), 
                               .SDcols = paste(substr(currency_names, 1, 3), rate, interest_name, sep = '_')]),
                  t(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)), 
                               .SDcols = paste(substr(currency_names, 5, 7), rate, interest_name, sep = '_')]),
                  t(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)),
                               .SDcols = paste(currency_names, rate, 'spot_log_basis', interest_name, sep = '_')]))
    temp <- cbind(temp, temp[, 1] - temp[, 2])
    rownames(temp) <- paste(currency_names, rate, interest_name, sep = '_')
    colnames(temp) <- c('funding_currency_interest', 'investing_currency_interest', 'avg_basis', 'difference')
    
    list_interest[[counter]] <- temp
    counter <- counter + 1
  }
}
names(list_interest) <- c('ois_1M', 'ois_3M', 'ibor_3M')

# Average interest rate, basis, correlations of basis and FX with S&P for single currency against USD.
# OIS 3M only.
# Weekly correlation taken on Wednesdays.  4 missing S&P due to 4th of July, Christmas, New Year.
sp_old <- fread(paste0(script_path, 'DataSkeleton/AssetClasses/SPX_dailypricehistory.csv'))
sp_new <- fread(paste0(script_path, 'DataSkeleton/AssetClasses/SPX_dailypricehistory_extended.csv'))
sp_old[, `:=` (Date = ymd(Date))]
sp_new[, `:=` (Date = ymd(DATE),
               SPX = as.numeric(SP500))]
sp <- rbind(sp_old, sp_new[Date >= ymd(max(sp_old$Date)), .(Date, SPX)])
rates_temp[, sp_level := sp$SPX[match(rates_temp$Date, sp$Date)]]

for (rate in c('ois')) {
  currency_names_single <- paste(single_table$fund_currency, single_table$invest_currency, sep = '_')
  interest_name_list <- c('3M')
  
  for (interest_name in interest_name_list) {
    rates_temp[, day_of_week := weekdays(Date)]
    for (c in 1:length(currency_names_single)) {
      rates_temp[, paste(currency_names_single[c], rate, 'log_1M_fwd_less_spot_basis', interest_name, sep = '_') := .SD[, 1] - .SD[, 2],
                 .SDcols = c(paste(currency_names_single[c], rate, 'log_basis_1M_fwd', interest_name, sep = '_'),
                             paste(currency_names_single[c], rate, 'spot_log_basis', interest_name, sep = '_'))]
    }
    rates_temp_filter <- rates_temp[day_of_week == 'Wednesday']
    rates_temp_filter[, paste(currency_names_single, rate, 'log_basis_change_3M', sep = '_') := lapply(.SD, function(x) x - shift(x, type = 'lag')),
                      .SDcols = paste(currency_names_single, rate, 'spot_log_basis_3M', sep = '_')]
    rates_temp_filter[, paste(substr(currency_names_single, 1, 3), 'log_fx_change_spot', sep = '_') := lapply(.SD, function(x) log(shift(x, type = 'lag')) - log(x)),
                      .SDcols = paste(substr(currency_names_single, 1, 3), 'fx_spot', sep = '_')] #Convention for FX return is long foreign short USD, so appreciation (lower tmr) is gain
    rates_temp_filter[, sp_return := log(sp_level) - log(shift(sp_level, type = 'lag'))]
    
    ## averages are calculated from rates_temp data table; correlations are from rates_temp_filter
    results_table <- data.table(pair = currency_names_single,
                                funding_currency_interest = unlist(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)), 
                                                                              .SDcols = paste(substr(currency_names_single, 1, 3), rate, interest_name, sep = '_')]),
                                invest_currency_interest = unlist(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)), 
                                                                             .SDcols = paste(substr(currency_names_single, 5, 7), rate, interest_name, sep = '_')]),
                                average_basis = unlist(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)),
                                                                  .SDcols = paste(currency_names_single, rate, 'spot_log_basis', interest_name, sep = '_')]),
                                sp_corr_fx = unlist(rates_temp_filter[, lapply(.SD, function(x) cor(x, sp_return, use = 'pairwise.complete.obs')),
                                                                      .SDcols = paste(substr(currency_names_single, 1, 3), 'log_fx_change_spot', sep = '_')]),
                                sp_corr_basis = unlist(rates_temp_filter[, lapply(.SD, function(x) cor(x, sp_return, use = 'pairwise.complete.obs')),
                                                                         .SDcols = paste(currency_names_single, rate, 'log_basis_change_3M', sep = '_')]),
                                average_return = unlist(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)),
                                                                   .SDcols = paste(currency_names_single, rate, 'log_return_1M_fwd', interest_name, sep = '_')]),
                                average_slope = unlist(rates_temp[, lapply(.SD, function(x) mean(x, na.rm = T)),
                                                                  .SDcols = paste(currency_names_single, rate, 'log_1M_fwd_less_spot_basis', interest_name, sep = '_')])
    )
    results_table[, interest_difference := funding_currency_interest - invest_currency_interest]
    corr_table <- results_table[, .(pair, average_return, average_slope, average_basis, interest_difference, sp_corr_fx, sp_corr_basis)]
    setorder(corr_table, -average_slope)
    
    rank_table <- data.table(pair = corr_table[, pair], 
                             average_return = frank(-corr_table$average_return),
                             average_slope = frank(-corr_table$average_slope),
                             average_basis = frank(-corr_table$average_basis), 
                             interest_difference = frank(-corr_table$interest_difference),
                             sp_corr_fx = frank(corr_table$sp_corr_fx), 
                             sp_corr_basis = frank(corr_table$sp_corr_basis))
  }
}

#---------------------------------------------------------
# Portfolio returns - 1M-forward 3M, 3M-forward 3M
# ** Note: annualized return & SR **
#---------------------------------------------------------

list_portfolio_return <- list()
list_portfolio_mean_star <- list()
list_portfolio_Sharpe_star <- list()

list_portfolio_return_covid <- list()
list_portfolio_mean_star_covid <- list()
list_portfolio_Sharpe_star_covid <- list()

list_portfolio_return_slr <- list()
list_portfolio_mean_star_slr <- list()
list_portfolio_Sharpe_star_slr <- list()

counter <- 1
for (rate in c('ois', 'ibor')) {
  
  if (rate == 'ois') {
    return_name_list <- c('1M_fwd_1M', '1M_fwd_3M', '3M_fwd_3M')
  } else if (rate == 'ibor') {
    return_name_list <- c('1M_fwd_3M', '3M_fwd_3M')
  }
  
  for (return_tenor in return_name_list) {
    portfolio_return_names <- names(basis_and_returns)[grepl(paste('portfolio', rate, 'log_return', sep = '_'), names(basis_and_returns)) &
                                                         grepl(return_tenor, names(basis_and_returns))]
    weighting_scheme_names <- sapply(portfolio_return_names, function(x) paste(strsplit(x, split = '_')[[1]][5], strsplit(x, split = '_')[[1]][6]))
    
    # store summary stats
    temp_matrix <- data.frame(matrix(, nrow = length(return_matrix_row) * length(portfolio_return_names), ncol = length(period_order)))
    colnames(temp_matrix) <- c(period_order)
    rownames(temp_matrix) <- paste(sapply(weighting_scheme_names, function(x) rep(x, length(return_matrix_row))), return_matrix_row)
    
    temp_matrix_covid <- data.frame(matrix(, nrow = length(return_matrix_row) * length(portfolio_return_names), ncol = length(covid_order)))
    colnames(temp_matrix_covid) <- c(covid_order)
    rownames(temp_matrix_covid) <- paste(sapply(weighting_scheme_names, function(x) rep(x, length(return_matrix_row))), return_matrix_row)
    
    temp_matrix_slr <- data.frame(matrix(, nrow = length(return_matrix_row) * length(portfolio_return_names), ncol = length(slr_order)))
    colnames(temp_matrix_slr) <- c(slr_order)
    rownames(temp_matrix_slr) <- paste(sapply(weighting_scheme_names, function(x) rep(x, length(return_matrix_row))), return_matrix_row)
    
    for (i in 1:length(portfolio_return_names)) {
      for (p in 1:length(period_order)) {
        temp_summary <- try(unlist(summary_stats_returns(temp_data_table = basis_and_returns[period == period_order[p]], column_to_summ = portfolio_return_names[i], 
                                                         mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                                         basis_or_percentage = 'basis')))
        if (class(temp_summary) == 'try-error') {
          temp_matrix[grepl(weighting_scheme_names[i], rownames(temp_matrix)), period_order[p]] <- NA
        } else {
          temp_matrix[grepl(weighting_scheme_names[i], rownames(temp_matrix)), period_order[p]] <- temp_summary
        }
      }
      for (q in 1:length(covid_order)) {
        temp_summary <- tryCatch(unlist(summary_stats_returns(temp_data_table = basis_and_returns[covid_period == covid_order[q]], column_to_summ = portfolio_return_names[i], 
                                                         mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                                         basis_or_percentage = 'basis')),
                                 error = function(e) {}) #error from no CHF post-2020
        if (is.null(temp_summary)) {
          temp_matrix_covid[grepl(weighting_scheme_names[i], rownames(temp_matrix_covid)), covid_order[q]] <- NA
        } else {
          temp_matrix_covid[grepl(weighting_scheme_names[i], rownames(temp_matrix_covid)), covid_order[q]] <- temp_summary
        }
      }
      for (q in 1:length(slr_order)) {
        temp_summary <- try(unlist(summary_stats_returns(temp_data_table = basis_and_returns[slr_period == slr_order[q]], column_to_summ = portfolio_return_names[i], 
                                                         mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                                         basis_or_percentage = 'basis')))
        if (class(temp_summary) == 'try-error') {
          temp_matrix_slr[grepl(weighting_scheme_names[i], rownames(temp_matrix_slr)), slr_order[q]] <- NA
        } else {
          temp_matrix_slr[grepl(weighting_scheme_names[i], rownames(temp_matrix_slr)), slr_order[q]] <- temp_summary
        }
      }
    }
    
    temp_title <- paste('Summary stats on returns from portfolios of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_portfolio_return[[counter]] <- temp_matrix
    
    temp_title <- paste('COVID stats on returns from portfolios of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_portfolio_return_covid[[counter]] <- temp_matrix_covid

    temp_title <- paste('SLR stats on returns from portfolios of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_portfolio_return_slr[[counter]] <- temp_matrix_slr
    
    # manually create stars
    star_mean <- matrix(, nrow = length(weighting_scheme_names), ncol = length(period_order))
    colnames(star_mean) <- period_order
    rownames(star_mean) <- weighting_scheme_names
    
    star_Sharpe <- matrix(, nrow = length(weighting_scheme_names), ncol = length(period_order))
    colnames(star_Sharpe) <- period_order
    rownames(star_Sharpe) <- weighting_scheme_names
    
    star_mean_covid <- matrix(, nrow = length(weighting_scheme_names), ncol = length(covid_order))
    colnames(star_mean_covid) <- covid_order
    rownames(star_mean_covid) <- weighting_scheme_names
    
    star_Sharpe_covid <- matrix(, nrow = length(weighting_scheme_names), ncol = length(covid_order))
    colnames(star_Sharpe_covid) <- covid_order
    rownames(star_Sharpe_covid) <- weighting_scheme_names
    
    star_mean_slr <- matrix(, nrow = length(weighting_scheme_names), ncol = length(slr_order))
    colnames(star_mean_slr) <- slr_order
    rownames(star_mean_slr) <- weighting_scheme_names
    
    star_Sharpe_slr <- matrix(, nrow = length(weighting_scheme_names), ncol = length(slr_order))
    colnames(star_Sharpe_slr) <- slr_order
    rownames(star_Sharpe_slr) <- weighting_scheme_names
    
    for (c in 1:length(weighting_scheme_names)) {
      star_mean[weighting_scheme_names[c], ] <- make_stars(values = list_portfolio_return[[counter]][paste(weighting_scheme_names[c], 't_mean'), ],
                                                           df = list_portfolio_return[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
      star_Sharpe[weighting_scheme_names[c], ] <- make_stars(values = list_portfolio_return[[counter]][paste(weighting_scheme_names[c], 't_Sharpe'), ],
                                                             df = list_portfolio_return[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
      star_mean_covid[weighting_scheme_names[c], ] <- make_stars(values = list_portfolio_return_covid[[counter]][paste(weighting_scheme_names[c], 't_mean'), ],
                                                           df = list_portfolio_return_covid[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
      star_Sharpe_covid[weighting_scheme_names[c], ] <- make_stars(values = list_portfolio_return_covid[[counter]][paste(weighting_scheme_names[c], 't_Sharpe'), ],
                                                             df = list_portfolio_return_covid[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
      star_mean_slr[weighting_scheme_names[c], ] <- make_stars(values = list_portfolio_return_slr[[counter]][paste(weighting_scheme_names[c], 't_mean'), ],
                                                                 df = list_portfolio_return_slr[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
      star_Sharpe_slr[weighting_scheme_names[c], ] <- make_stars(values = list_portfolio_return_slr[[counter]][paste(weighting_scheme_names[c], 't_Sharpe'), ],
                                                                   df = list_portfolio_return_slr[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
    }
    
    list_portfolio_mean_star[[counter]] <- star_mean
    list_portfolio_Sharpe_star[[counter]] <- star_Sharpe
    list_portfolio_mean_star_covid[[counter]] <- star_mean_covid
    list_portfolio_Sharpe_star_covid[[counter]] <- star_Sharpe_covid
    list_portfolio_mean_star_slr[[counter]] <- star_mean_slr
    list_portfolio_Sharpe_star_slr[[counter]] <- star_Sharpe_slr
    
    counter <- counter + 1
  }
}
names(list_portfolio_return) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_portfolio_mean_star) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_portfolio_Sharpe_star) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

names(list_portfolio_return_covid) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_portfolio_mean_star_covid) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_portfolio_Sharpe_star_covid) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

names(list_portfolio_return_slr) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_portfolio_mean_star_slr) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_portfolio_Sharpe_star_slr) <- c('ois_1M_fwd_1M', apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

#---------------------------------------------------------
# Monthly portfolio returns - 1M-forward 3M, 3M-forward 3M
# ** Note: annualized return & SR **
#---------------------------------------------------------

## Construct monthly returns
portfolio_names <- c('classic_carry', 'single_3v3', 'equal_dollar', 'equal_pairs', 'top5_basis', 'top6_pc')
list_temp <- list()
counter <- 1
for (rate in c('ois', 'ibor')) {
  for (return_tenor in c('1M_fwd_3M', '3M_fwd_3M')) {
    for (port_name in portfolio_names) {
      temp <- basis_and_returns[, c('Date', paste('portfolio', rate, 'log_return', port_name, return_tenor, sep = '_')), with = FALSE] 
      temp <- temp[complete.cases(temp)]
      temp[, paste('dd_port', rate, substr(port_name, 1, regexpr('_', port_name) - 1), return_tenor, sep = '_') := max(Date), by = floor_date(Date, unit = 'month')]
      temp <- temp[Date <= trading_end_date]
      temp <- temp[Date == get(paste('dd_port', rate, substr(port_name, 1, regexpr('_', port_name) - 1), return_tenor, sep = '_'))]
      temp <- temp[, Date := as.integer(paste0(year(get(paste('dd_port', rate, substr(port_name, 1, regexpr('_', port_name) - 1), return_tenor, sep = '_'))), 
                                               sprintf('%02d', month(get(paste('dd_port', rate, substr(port_name, 1, regexpr('_', port_name) - 1), return_tenor, sep = '_'))))))]
      list_temp[[counter]] <- temp
      counter <- counter + 1
    } 
  }
}
return_monthly <- Reduce(function(x, y) {merge(x, y, all = TRUE, by = 'Date')}, list_temp)
return_monthly[, period := ifelse(Date < 200707, 'Pre-Crisis',
                                  ifelse(Date < 201007, 'Crisis', 'Post-Crisis'))]

## Summarize
list_monthly_return <- list()
list_monthly_mean_star <- list()
list_monthly_Sharpe_star <- list()
counter <- 1
for (rate in c('ois', 'ibor')) {
  return_name_list <- c('1M_fwd_3M', '3M_fwd_3M')
  
  for (return_tenor in return_name_list) {
    portfolio_return_names <- names(return_monthly)[grepl(paste('portfolio', rate, 'log_return', sep = '_'), names(return_monthly)) &
                                                      grepl(return_tenor, names(return_monthly))]
    weighting_scheme_names <- sapply(portfolio_return_names, function(x) paste(strsplit(x, split = '_')[[1]][5], strsplit(x, split = '_')[[1]][6]))
    
    # store summary stats
    temp_matrix <- data.frame(matrix(, nrow = length(return_matrix_row) * length(portfolio_return_names), ncol = length(period_order)))
    colnames(temp_matrix) <- c(period_order)
    rownames(temp_matrix) <- paste(sapply(weighting_scheme_names, function(x) rep(x, length(return_matrix_row))), return_matrix_row)
    
    for (i in 1:length(portfolio_return_names)) {
      for (p in 1:length(period_order)) {
        temp_matrix[grepl(weighting_scheme_names[i], rownames(temp_matrix)), period_order[p]] <- 
          unlist(summary_stats_returns(temp_data_table = return_monthly[period == period_order[p]], column_to_summ = portfolio_return_names[i], 
                                       mths_in_return = as.numeric(substr(return_tenor, 0, 1)), mths_in_contract = as.numeric(substr(return_tenor, 8, 8)),
                                       basis_or_percentage = 'basis', NeweyWest = F)) #NON-overlapping returns
      }
    }
    
    temp_title <- paste('Summary stats on returns from MONTHLY portfolios of', rate, unlist(strsplit(return_tenor, split = '_'))[[1]], 'fwd', unlist(strsplit(return_tenor, split = '_'))[[3]])
    list_monthly_return[[counter]] <- temp_matrix
    
    # manually create stars
    star_mean <- matrix(, nrow = length(weighting_scheme_names), ncol = length(period_order))
    colnames(star_mean) <- period_order
    rownames(star_mean) <- weighting_scheme_names
    
    star_Sharpe <- matrix(, nrow = length(weighting_scheme_names), ncol = length(period_order))
    colnames(star_Sharpe) <- period_order
    rownames(star_Sharpe) <- weighting_scheme_names
    
    for (c in 1:length(weighting_scheme_names)) {
      star_mean[weighting_scheme_names[c], ] <- make_stars(values = list_monthly_return[[counter]][paste(weighting_scheme_names[c], 't_mean'), ], 
                                                           df = list_monthly_return[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
      star_Sharpe[weighting_scheme_names[c], ] <- make_stars(values = list_monthly_return[[counter]][paste(weighting_scheme_names[c], 't_Sharpe'), ],
                                                             df = list_monthly_return[[counter]][paste(weighting_scheme_names[c], 'sample_size'), ])
    }
    
    list_monthly_mean_star[[counter]] <- star_mean
    list_monthly_Sharpe_star[[counter]] <- star_Sharpe
    
    counter <- counter + 1
  }
}
names(list_monthly_return) <- c(apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_monthly_mean_star) <- c(apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))
names(list_monthly_Sharpe_star) <- c(apply(expand.grid(c('1M_fwd_3M', '3M_fwd_3M'), c('ois', 'ibor')), 1, function(tuple) paste(tuple[2], tuple[1], sep = '_')))

## Save
save(list_single_return, list_single_mean_star, list_single_Sharpe_star, list_interest, corr_table, rank_table,
     list_portfolio_return, list_portfolio_mean_star, list_portfolio_Sharpe_star,
     list_monthly_return, list_monthly_mean_star, list_monthly_Sharpe_star,
     list_single_return_covid, list_single_mean_star_covid, list_single_Sharpe_star_covid, 
     list_portfolio_return_covid, list_portfolio_mean_star_covid, list_portfolio_Sharpe_star_covid,
     list_single_return_slr, list_single_mean_star_slr, list_single_Sharpe_star_slr, 
     list_portfolio_return_slr, list_portfolio_mean_star_slr, list_portfolio_Sharpe_star_slr,
     file = paste0(script_path, 'OutputInterim/return_summary.RData'))

